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Abstract - We present a new Lattice Boltzmann (LB) formulation to solve the Maxwell equations 
for electromagnetic (EM) waves propagating in a heterogeneous medium. By using a pseudo- vector 
discrete Boltzmann distribution, the scheme is shown to reproduce the continuum Maxwell equa- 
tions. The technique compares well with a pseudo-spectral method at solving for two-dimensional 
wave propagation in a heterogeneous medium, which by design contains substantial contrasts in 
the refractive index. The extension to three dimensions follows naturally and, owing to the recog- 
nized efficiency of LB schemes for parallel computation in irregular geometries, it gives a powerful 
method to numerically simulate a wide range of problems involving EM wave propagation in 
complex media. 
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Introduction. — Over the last two decades, the LB 
method in the Bhatnagar-Gross-Krook (BGK) approxi- 
mation (see pQ) has met with significant success in simu- 
lating a broad spectrum of complex flow phenomena [2,3, 
ranging from low-Reynolds flows in porous media to fully 
developed turbulent flows in complex geometries, multi- 
phase flows and, more recently, relativistic flows [3]. The 
application of LB techniques to EM wave propagation phe- 
nomena is comparatively far less developed. Given the 
paramount role of EM phenomena in science and technol- 
ogy, including fast-emerging applications such as metama- 
terials and transformation optics [S] , it is of great interest 
to investigate whether LB techniques are able to improve 
simulations of complex EM phenomena as they have for 
fluid flows. 

LB schemes for wave propagation have been proposed 
before in literature, [SHE], but it is only recently that such 
schemes have been specifically tailored to EM phenomena. 
To this end, a basic issue had to be addressed, namely the 
fact that, unlike hydrodynamics, EM interactions are gov- 
erned by an anta-symmetric field tensor, a structure that 
does not naturally emerge from standard kinetic theory. 
To circumvent this problem, Mendoza et al. [IIHHI] pro- 
posed a scheme wherein the electric and magnetic fields are 



represented by separate discrete Boltzmann distributions, 
each moving along distinct lattice directions. Although 
an important conceptual advance, the resulting scheme 
appears computationally intensive. 

A more straightforward formulation for EM wave prop- 
agation in plasmas has been recently proposed by Dol- 
lar ( [T2], [13] )• By promoting the discrete Boltzmann 
distribution from a scalar to a vector quantity (also see 
[7], [9]), and expressing the curl of electric field in di- 
vergence form, Dellar developed an elegant and compact 
scheme in which magnetic field emerges as the "vector 
density" associated with the vector-valued discrete Boltz- 
mann distribution. The Maxwell equations emerge when 
a Hermite-projection procedure is applied to the resulting 
kinetic equation. However, antisymmetry of the relevant 
EM tensor is not preserved in time, so it must be enforced 
at each time-step. 

Motivated by Dollar's core idea of the discrete Boltz- 
mann distribution not having to be a scalar, we formu- 
late a new scheme that is computationally inexpensive, 
inherently maintaining antisymmetry of the EM tensor. 
To this end, we introduce a tensorial distribution func- 
tion, with built-in antisymmetry g a $ = —gp a , where greek 
subscripts a, /? run over spatial dimensions. Tensor g a p 
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is in fact a pseudo-vector, entailing only D independent 
components in _D-dimensional space. Furthermore, and 
most importantly for practical applications, we develop a 
new formulation wherein heterogeneity is realized through 
space-time-depcndcnt permittivities and light speed, and 
embedded within a source term in the Maxwell equations. 
The scheme developed here permits us to address the im- 
portant problem of EM wave propagation in strongly het- 
erogeneous materials. The scheme is validated by direct 
comparison with a pseudo-spectral method for the case of 
wave propagation in both homogeneous and heterogeneous 
media. 

The Maxwell Equations. — We start with the 
Maxwell equations defined in arbitrarily heterogeneous, 
charge-free media. In principle, we allow for frequency- 
dependent dielectric permittivity, which introduces a tem- 
poral convolution, resulting in integro-differential Maxwell 
equations. For the sake of simplicity, we consider magnetic 
permeability as only a function of space (and therefore 
constant in time); this assumption may be relaxed and 
frequency-dependent permeability may be treated with- 
out further conceptual difficulty. The Maxwell equations 
read as follows: 

dtB a = —e a p 1 dpE 1 (1) 
dtE a = t a p 1 c 2 dpB 1 + S a , (2) 
S a = c 2 e a j3y Bp d 7 In // 

- J —- ( dt'e{*,t-t')EJ*,t'), (3) 
eo Jo 

d a B a = 0, (4) 

where wc use Einstein's summation convention, e a ^ 7 is the 
Levi-Civita tensor, E a is the electric field, B a the mag- 
netic field, j a the current, /i (x), e ( x ), e( x , 0: are P er ~ 
mittivities that may be inhomogeneous, e = d t e, x and t 
are spatial and temporal coordinates respectively, speed 
of light c(x) = 1/ y/JIoeo, and S a (x, t) is the source term 
into which all inhomogeneities are collected. The temporal 
convolution term is a causal model of frequency-dispersive 
permittivity and is termed the polarization field. In deriv- 
ing these equations, we have split the electric displacement 
field, which is a sum of electric and polarization fields, into 
two terms and treat the latter as a source. We do so in 
order to bring the Maxwell equations into a conservative 
form, amenable to solution by LB methods. 

Various numerical techniques may be utilized to solve 
these integro-differential equations - one of the most popu- 
lar is to treat the convolution in terms of auxiliary differen- 
tial equations (ADEs) (e.g., [E]). As in standard control 
theory, the Laplace transform of e(x, t) may be written in 
terms of its zeros and poles [15] and using these, a system 
of differential equations, whose effective continuous-time 
response is the convolution, may be constructed. A sim- 
ilar method may be used in order to address frequency- 
dispersive magnetic permeabilities. 



Although conceptually straightforward, the computa- 
tional treatment of these convolutions is rather laborious, 
and consequently, we shall leave it for a future separate 
study. In what follows, we focus on a simpler case in 
which dielectric and magnetic permittivities are allowed 
to vary spatially and but are temporally stationary. 

A key step to developing an LB formulation is to cast 
the Maxwell equations in conservative form. To this end, 
we express the curl of electric field as the divergence of 
an anti-symmetric second-order tensor, h a p = — e 7Q , J g£' 7 , 
such that dpAp a = e a 0jdpEj. As a result, the induction 
equation (TJJ) reduces to d t B a + dpKp a = 0. By invoking 
the Levi-Civita identity, \e ia p e™^ = S^ T , the electric 
field may be recovered through: E 1 = — n^yafi^-afi- In 
terms of new variables B a and A a p , the Maxwell equations 
take on the following conservative form: 

d t B a +d p A Pa = 0, (5) 
d t A al3 + c 2 n a p = -e 7Q(3 S' 7 , (6) 

where Q a p = d a Bp — dpB a is the curl of the magnetic 
field (the wedge curl V A B in Clifford algebra termi- 
nology). The wedge curl of the magnetic field may also 
be expressed as the divergence of a triple-rank tensor, so 
that the Maxwell equations appear in fully conservative 
form but we are able to avoid this additional complica- 
tion. Note that if d a B a = at t = 0, then it will re- 
main so for all time since the partial derivative with re- 
spect to a of equation ([5]) gives d t d a B a = —d a dpKp a , or, 
d t d a B a = —e^padadpEj = 0. 

Lattice BGK formulation. The next step is to 
formulate a discrete kinetic equation, whose continuum 
limit reproduces the Maxwell equations in the form (O, © 
given above. To this end, we define the pseudo-vector 
distribution function gi = gi a /3, where (a,/3) are vector 
indices and i labels the discrete identity of the particle 
moving with velocity Cj, to be detailed shortly. 

The pseudo- vector distribution is postulated to obey the 
LB equation in BGK form 

A lgl = g i (x + c i At,t + At) -ft(x,f) 

- (°) 

= -— — — — At + T,-Ai, (7) 

T 

where we define the tensor source term as: Tj = Ti a p = 
— ^c-yapSy. The local equilibrium is chosen so as to re- 
produce equations ©, (O in the continuum, and it reads 
as follows: 

9i% = I A "/J + c "±Bp - CipB a ] , (8) 

where ci a and Wi arc discrete particle velocities and lat- 
tice weights, respectively. The structural difference with 
hydrodynamics is apparent here: inner scalar products are 
replaced by outer (wedge) products. 

In D-dimcnsional space, these equations are formulated 
in a nearest-neighbor lattice with 2D discrete speeds, plus 
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a rest particle. For instance, in two spatial dimensions 
and counterclockwise counting, we have ci = (1,0), C2 = 
(0,1), c 3 = (-1,0), c 4 = (0,-1) and c = (0,0). The non- 
moving rest particle '0' (co = (0,0)), is instrumental in 
implementing a spatially varying wave propagation speed. 
To this end, weights of moving particles are all set to the 
value i«i(x) = (1 — wq(x))/2D. 

According to standard LB theory (e.g., [7J, [5]), wave 
propagation speed (squared) is given by the weighted sum 
of squared particle speeds, c 2 = J2i w i c i- As a rcsu lt, for 
the 6 + 1 speeds lattice, we obtain c 2 (x) = (1 — wq(x.))/3, 
which is the desired spatial dependence of lattice light 
speed. 

Note that by choosing different weights for particles 
moving along different directions x,y and z, we are able to 
model anisotropic media. 

The first three moments of the equilibrium distribution 
function (Eq. [8]) are readily computed 

i 

G L°/3 7 = Yl *l9?2p = S ~ta B P - 5 lP B a, (10) 
i 

G S 7 « = Yl c ^ c ^& = <vw (11) 

i 

We also stipulate that the distribution function (Eq. [5]) 
satisfies "mass-momentum" conservation constraints 

EGfotf-ffiXte = 0, (12) 

i 

with fa = (l,Cj 7 ). 

The above constraint implies that the non-equilibrum 
component of gi a p must contribute zero change to local 
mass and momentum, which is a distinctive feature of 
model BGK equations in general. 

Higher-order conservations may also be enforced, but 
these require the use of more complex lattices, with higher 
symmetries. Although certainly within the realm of pos- 
sibility, it is more labor intensive, especially in connection 
with complex geometries. 

To analyze the continuum limit, we Taylor expand the 
left side (streaming operator) of equation ([7]) to first order 
in At, to obtain: A, ~ AtDi = At(dt + Cj • V), where 
Di is the Lagrangian derivative along the i-th discrete di- 
rection. Summing over lattice speeds and invoking con- 
straints (HI, 

dtG a j3 + <9 7 G Q( 3 7 = — -^e ja ^S^, (13) 
dtGapj + d K G a p 1K = 0, (14) 

where, by definition, G a p = J2i9iap, G a p~ ( = J2i9ia^c i7 , 
G a p 1K = J2i9iai3C ly c iK . By construction, G a p = G^°j 
and G Q( 3 7 = G^l , which, along with the approximation 



G a p 7K ~ G a p~ K i and accounting for ©, (|TU1) . pT|) . allows 
us to rewrite moment equations (|13[). p4H 

d t A a p + c 2 (d a Bp - dpB a ) = -e 7Q/3 S' 7 , (15) 
S^adiBp - 5 K j3d t B a + S^djAap = 0. (16) 

It may be verified that for /3 = k ^ a (the case 
a = ft being trivial) returns (O. In summary, moment 
equations (|15[) . ([TB"]) are shown to reduce to the Maxwell 
equations in conservative form. 

As is well known in the case of fluids, a consistent analy- 
sis of dissipative terms requires the streaming operator to 
be expanded to second order, i.e., Di = di + ^di. Lengthy 
algebra shows that the right side of equation ([16]) acquires 
a dissipative term of the form c 2 (r — At/2)A(d/3A a p). As 
discussed in [7J for the case of scalar waves, this is set to 
zero by choosing r = At/ 2 (r = 1/2 in lattice units). 

The above formalism readily translates into a simple 
and actionable algorithm, consisting of two basic steps: 
collision and streaming. In the former, one first prepares 
the so-called post-collisional state as follows: 

gLpi*' t) = C 1 - — )3ia/3(x, t) + —g-L(^ *) 

T T r 

+ T iaP ( x ,t)At, (17) 

where macroscopic variables such as electric and magnetic 
fields, needed to compute local equilibrium, are obtained 
from moment equations (O and (|10j) . Subsequently, in the 
streaming step, the post-collisional state is simply shifted 
to a neighboring lattice point depending on the direction 
and sign of the velocity, namely: 

9ia p{x + Ci At,t + At) = g' iafs (x,t) (18) 

It may be verified that equations (|17p and (|T5|) are strictly 
equivalent to the LB equation ©. 

The resulting numerical procedure is elegant and easy 
to code, as thoroughly discussed in numerous introductory 
texts on this topic (e.g., [2"0]). 

Numerical results: 2D wave propagation. 

Since one of the highlights of our scheme is built-in 
antisymmetry, we first inspect numerical errors in main- 
taining a divergence-free magnetic field in a homogeneous 
medium, where c 2 (x,y) = 1, the lattice speed being mea- 
sured in units of its uniform value. Particle streaming and 
collisions are performed for all components of the distri- 
bution function on a 4 + 1-speed two-dimensional lattice, 
and we operate at r = 1/2 (in lattice units A< ~ 1). The 
implementation is akin to standard LB (e.g., [3]), the only 
difference being the inclusion of a multi-component distri- 
bution function. 

Homogeneous media. We excite waves by forc- 
ing the medium through current density j z (x,y,t) = 
exp [-[(x - 0.3) 2 + (y- 0.4) 2 ]/(2 x 0.03 2 )] , where e = 1, 
and (x,y) £ [0,1). All fields are initialized to zero. The 
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general wavenumber (k 
described by 



\k Xl ky\) dependent error e is 



Jdk' J"(k,k') \ik' -B(k') 
Jdk' T(k,k') |B(k')| 



(19) 



where we set J- = 1 to obtain error in the L2 norm, 
T = exp (— (|k| — |k'|) 2 /2fj 2 ) to model the variation of 
error as a function of absolute wavenumber, and thirdly, 
T = cxp (— |k — k'| 2 /2er 2 ) to capture error dependence 
on angle of propagation at fixed absolute wavenumber. 
These forms of error are plotted in Figure [TJ In gen- 
eral, we expect measures of model accuracy to reflect a 
second-order convergence rate. For example, the error in 
enforcing Gauss' law of charge conservation is also likely 
to be accurate only to second order (although not con- 
firmed by these tests). In the 2D cases considered here, 
the divergence of the electric field is identically zero here 
since E = {0, 0, E z {x, y)} -> VE = d z E z = 0. 

Heterogeneous media. Next, we simulate wave 
propagation in an inhomogeneous medium, with 
the initial condition B x = (y — 0.35)/(x, y), 
B y = -(x - 0.5)f(x,y), where f(x,y) = 
lOexp {-[(x - 0.5) 2 + (y - 0.35) 2 ]/(2 x 0.04 2 )}. As 
may be seen on the lower-left panel of Figure the 
shape, size and magnitude of dielectric "defects" in this 
calculation bear a qualitative resemblance with the inter- 
mediate matched-impedance zero-index material region 
(MIZIM) discussed in [TBJ. However, different from [15] , 
who consider frequency-dependent permittivity (implying 
a convolution with electric field), we only consider the 
time-stationary case. Also, we do not include inlet and 
outlet vacuum regions, but implement periodic boundary 
conditions. In order to avoid finite-size effects due to 
recirculating waves, the simulation is terminated before 
waves reach boundaries. 

The initial Gaussian wave-packet splits into up- and 
downward propagating components. Waves of finite spa- 
tial extent refract into (away from) regions of low (high) c, 
because the portion of the waveform within the inhomo- 
geneity propagates comparatively slower (faster). More- 
over, since we study linear waves, wave frequency may be 
regarded as invariant, so that wavelength A oc c. Conse- 
quently, waves exhibit locally smaller (larger) wavelengths 
in regions of low (high) c. This implies that energy, de- 
fined as £ = J2 a (^a + E^/c 2 ), tends to concentrate in 
regions of low c, because waves spend larger fractions of 
time in these areas. The (x, y) components of the Poynt- 
ing energy-flux vector, V a — e a p 1 EpB 1 , and electric field 
E z are also shown. These interpretations arc confirmed by 
visual inspection of properties of the wavefield and c 2 (x) , 
as shown in Figures [2] and |3l Note that the magnitude of 
c 2 in vacuum is necessarily greater than the largest value 
in this medium. The refractive index, n cx c _1 , is seen to 
vary by a factor of four (c 2 ~ 0.1 — 1.6), comparable to 
contrasts used in [1^ to fine-tune transmission (reflection) 
coefficients across the MIZIM region. 
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Fig. 1: Errors in maintaining a divergenceless magnetic field as 
denned by equation (|19p . On the upper left panel is the error 
in the L2 norm plotted as a function of time; note that there 
is an initial transient in the magnetic field, due to the applied 
current, followed by a period where the error undergoes some 
fluctuations. On the upper-right panel, the error as a func- 
tion of grid spacing in the £2 norm is plotted, along with a 
nominal second-order convergence rate curve. For a given cal- 
culation (n x — 32 here), the error is Fourier-decomposed and 
plotted as a function of wavenumber or equivalently, points 
per wavelength on the lower-left panel. The nominal line is 
also plotted, but the error is seen to only approximately fall at 
a second-order rate. Lastly, the error is anisotropic, as shown 
on the lower-right panel. Distance from the coordinate center 
denotes error magnitude; waves propagating at 22.5° with re- 
spect to the axes are seen to be more inaccurately resolved than 
in other directions. This is a function of the choice of lattice 
and, in principle, error anisotropy may be reduced through the 
use of higher-order lattices. 



In order to test the accuracy of the LB solution, we 
repeat this calculation using a pseudo-spectral solver, in 
which spatial derivatives are computed spectrally and time 
stepping is achieved through the application of an opti- 
mized Rungc-Kutta scheme [T5]. In Figure we also 
show a comparison between outputs of the two simula- 
tions. Grids in both solvers contain 512 x 512 points; 
LB and pseudo-spectral solutions are very similar, with 
a 0.6% Li norm of the difference. The LB calculation 
at this resolution is approximately 4 times faster than its 
pseudo-spectral counterpart; however, these codes were 
not written with optimization and efficiency in mind and 
the number may only be interpreted loosely. 

Outlook. — In summary, we demonstrate that top- 
down-prescribed distribution functions of "particles" , not 
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Fig. 2: Upper panels show snapshots of B x and B y at t = 1.272, 
with the initial condition described in the text; the lower left 
shows the spatial distribution of c 2 . Axes of the three contour 
plots are x and y. In both calculations, the grid contained 
512 x 512 points and the domains were horizontally periodic, 
with x, y £ [0, 1). The horizontal colorbar applies to the upper 
panels while the vertical bar describes the range in c 2 . Regions 
of low wave-speed (compared with reference c 2 = 1; note that 
c 2 > 1 in vacuum, although we do not specify a value, since 
vacuum regions are not included in the computational domain) 
cause waves to refract towards them and vice-versa. The lower- 
right panel compares pseudo-spectral and LB solutions along 
a cut at constant x (location indicated by the dashed line in 
the B x plot). The technique is seen to accurately simulate 
wave propagation through regions where refractive index n (oc 
c _1 ) shows substantial local variations. The L2 norm of the 
difference between LB and pseudo-spectral solutions is 0.6%. 
In order to prevent boundary periodicity from affecting the 
results, a larger computational domain may be chosen and the 
calculation terminated before waves reach the boundaries. 



based on any known microscopic kinetic theory of EM, suc- 
ceed in reproducing the behavior predicted by the Maxwell 
equations. Moreover, our result provides a further exam- 
ple of lattice kinetic theory as an efficient tool to simulate 
continuum-physics phenomena via a particle-like formal- 
ism, well suited to handling complex geometries and show- 
ing excellent scalability on modern parallel computers [18] . 
The present method may be extended in many directions, 
such as invisibility cloaks, i.e., systems wherein the use 
of meta-materials allows for selected regions of space to 
become inaccessible to light, analogous to metric holes in 
space-time (e.g., Fig. 3 in [17]). Another possibility is to 
extend it to study plasma phenomena in confined geome- 
tries. 

We note potential limitations of this technique and de- 
scribe directions for future work. Like most LB meth- 
ods, the present scheme is formulated in a space-time uni- 
form lattice, so that space and time resolution may not 
be changed independently unless locally-adaptive formu- 
lations of the method are put in place. Such adaptive LB 
formulations do indeed exist [21] , although they are usu- 



Fig. 3: Upper panels show the x — y components of Poynting 
vector, and lower panels display energy and electric field E z . 
Axes of all plots are x and y, on a periodic grid of 512 x 512 
points and x, y £ [0, 1). Since waves spend more time in regions 
of low wave-speed (shorter wavelengths), energy appears to be 
concentrated in these areas. The top-right scale applies to both 
upper panels. 



ally significantly more laborious than the native version on 
uniform lattices. A second point regards numerical stabil- 
ity in the presence of sharp interfaces, resolved by just a 
few lattice points. At such sharp interfaces, higher order 
terms may no longer be negligible, thereby introducing un- 
wanted dissipativc effects. Efficient implementations with 
non-local permittivities, as mentioned earlier on in this pa- 
per, are likely to require a careful analysis of the auxiliary 
equations to be coupled with the native LB equation for 
waves. This is similar to the way auxiliary equations are 
coupled to the LB equation for the modeling of turbulent 
fluid flows [22] . 

The present LB scheme may be regarded as a special 
type of finite-difference method, in which streaming is ex- 
act and local conservations are built-in and accurate to 
machine round-off. Although the technique is just second- 
order in space-time, the above properties make the error 
prcfactors sufficiently small so as to render its performance 
competitive with higher-order methods, including spectral 
ones. 

From a mathematical stand point, the inclusion of 
source terms, reflecting external sources and/or inho- 
mogeneities, is straightforward, once the expression of 
these sources in the continuum is known. However, the 
strengths of such terms may place stringent constraints 
on the stability of the scheme, aspects that remain to be 
investigated. 

Although we have only considered periodic boundaries, 
we note that there exists literature on the implementation 
of other boundary conditions, such as reflecting, absorb- 
ing etc. (e.g., [3] and references therein) As a result, the 
formulation of different types of boundary conditions for 
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the wave-LB scheme might follow from previous develop- 
ments although such a possibility remains to be studied in 
full detail. 

Despite their importance, the above developments that 
do not challenge the basic merits of the scheme discussed 
in this paper. As a result, we believe that the LB scheme 
presented in this work should offer a fast and flexible com- 
putational tool to assist, complement and possibly even 
anticipate experimental research on invisibility cloaks and 
related phenomena in modern optics and photonics re- 
search |%P%]. 

* * * 
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